% Scripts written to read electric field data from
% Kakioka observatory, concatenate files and make 10 day
% time series of X and Y
% Dec 2013 (SFO)

clear

stations = ['kak';'mmb';'kny'];
%stations = ['kny'];

% First make observatory magnetic database and then get the electric field
% data for that period. 

load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_NovSq_2003.mat','fday');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Apr_2000.mat','fday');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Jul_2000.mat','fday');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Oct_2003.mat','fday');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Aug_2000.mat','fday');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Mar_2001.mat','fday');

    fday_start = fday(1) - (30/ ( 24 * 3600)); % adjust 30 sec
    fday_end = fday(end) - (30/ ( 24 * 3600));

    
for j = 1:size(stations,1),
st = 1;
en = 1440;
station = stations(j,:);


for i = floor(fday(1)) : floor(fday(end)),
    
    % daily minute e-field data from observatories
    fname = ['/Users/manojnair/data/obs_elec_data/' station '/' station datestr(i,'yyyymmdd') 'dgef.min'];

    
    data = read_iaga2002_1min_electric(fname);
    
    E_X_n(st:en,:) = data.MINUTELY(:,2);
    E_Y_n(st:en,:) = data.MINUTELY(:,3);
    fday_n(st:en) = data.MINUTELY(:,1);
    
    st = st + 1440;
    en = en + 1440;
    
end;
   

% get the data vector for the specific time
% + 10/(24*3600) to eleminates error from small round offs

L = fday_n  >= fday_start - 10/(24*3600) & fday_n <= ( fday_end + 10/(24*3600)) ;
E_X_n = E_X_n(L);
E_Y_n = E_Y_n(L);
fday_n = fday_n(L);


% correct the baseline jump for data KNY during Jul 2000

if (strcmp(datestr(fday(1),'yyyymmm'), '2001Mar') & strcmp(station,'kny'))
    
E_Y_n(10055:end) =  E_Y_n(10055:end) - (170-147.9);
E_X_n(10055:end) =  E_X_n(10055:end) - ( 231.4 - 209.3);

E_Y_n(10055:10261) = E_Y_n(10054);
E_X_n(10055:10261) = E_X_n(10054);

    
end;


% QC

L = isnan(E_X_n);
E_X_n = interp1(fday_n(~L), E_X_n(~L), fday_n,'linear',median(E_X_n(~L)));
L = isnan(E_Y_n);
E_Y_n = interp1(fday_n(~L), E_Y_n(~L), fday_n,'linear',median(E_X_n(~L)));

% 
% subplot(211);
% plot(E_X_n);
% subplot(212);
% plot (E_Y_n);
% title(station);
% 
% pause

% detrend

E_X_n = detrend(E_X_n);
E_Y_n = detrend(E_Y_n);




% add metadata

phi_gg(j) = data.LON;
theta_gg(j) = 90 - data.LAT;
names(j,:) = data.CODE;
E_X(:,j) = E_X_n';
E_Y(:,j) = E_Y_n';
fday = fday_n;

clearvars E_X_n E_Y_n st en fday_n L fname data 

end;

clearvars fday_start fday_end stations station i j  data
eval(['save /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Sq' datestr(fday(1),'mmmyy')]);

%% convert the Magnetic and electric field data for Nov 2003 from Alexei to common format

% Electric field data
load obsdata_elec.mat;
fday_c = fday + datenum(2000,1,1);
% center data on storm date
L = fday_c >= datenum(2003,11,17) & fday_c < datenum(2003,11,27);
E_X(~L,:) = [];
fday_c(~L) = [];
fday = fday_c;
E_Y(~L,:) = [];
names = name;
E_X = detrend(E_X);
E_Y = detrend(E_Y);
clearvars L name fday_c
eval(['save /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_' datestr(fday(1),'mmmyy')]);

clear;
% Magnetic data
load /Users/manojnair/projects/GIC/obs_nov_2003.mat
fday_c = fday + datenum(2000,1,1);
% center data on storm date
L = fday_c >= datenum(2003,11,17) & fday_c < datenum(2003,11,27);
X(~L,:) = [];
Y(~L,:) = [];
Z(~L,:) = [];
fday = fday_c(L);
clearvars L  fday_c

save /Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Nov_2003.mat




%% merge OBEM and Observatory data for 2000 Apr

load('/Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Apr00.mat')

E_X_n = E_X;
E_Y_n = E_Y;
fday_n = fday;
names_n = names;
phi_gg_n = phi_gg;
theta_gg_n = theta_gg;

load('/Users/manojnair/data/obs_elec_data/storm_selected_data/obsem_elec_apr00.mat')

theta_gg(3) = []; %OBEM 3
phi_gg(3) = [];

theta_gg = 90 - theta_gg;

E_X = [E_X' E_X_n];
E_Y = [E_Y' E_Y_n];
phi_gg = [phi_gg phi_gg_n];
theta_gg = [theta_gg theta_gg_n];
names = [names; names_n];

save /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Apr00.mat phi_gg theta_gg names E_X E_Y fday


% save the data



save  /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Apr00.mat E_X E_Y fday theta_gg phi_gg names



%% Testing the data sets

clear
close all

%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Apr_2000.mat');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Jul_2000.mat');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Oct_2003.mat');
%load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Aug_2000.mat');
load('/Users/manojnair/data/obs_mag_data/storm_selected_data/obs_Mar_2001.mat');

subplot(211);

plot(fday, nanmean(X),'r');
hold on
subplot(212);

plot(fday, nanmean(Y),'r');
hold on
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Apr00.mat  
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Jul00.mat  
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Oct03.mat  
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Aug00.mat  
%load /Users/manojnair/data/obs_elec_data/storm_selected_data/obs_elec_Mar01.mat  



subplot(211);

plot(fday, E_X);
hold on

subplot(212);

plot(fday, E_Y);
hold on



